##### Analytical results     #----------------------------------------------#

# This script generates all figures and regression tables that involve additional
# data sources such as V-Dem or the Oxford Covid-Tracker

#------------- Load and merge V-Dem data --------------------------------------------

# Load topic scores
load("output/acled_btm.Rda")

# Load V-Dem data (subset)
vdem_sub <- rio::import("data/vdem_subset.rda")

acled_topics_vdem <- left_join(acled_btm, vdem_sub, by = c("cowcode", "year")) %>% 
  select(data_id:corona_topic, max_topic, v2x_polyarchy:e_regionpol_6C, everything()) 


# Aggregtate at region level
acled_regions <- acled_topics_vdem %>% 
  group_by(e_regionpol_6C) %>% 
  summarise(across(c("health_care",
                     "restrictions_masks",
                     "vaccination",
                     "imprisonment_crime",
                     "mishandling_corruption",
                     "economic_consequences",
                     "business_restrictions",
                     "education",
                     "v2x_libdem"), 
                   ~ mean(., na.rm = T))) %>% 
  mutate(regio = case_when(e_regionpol_6C == 1 ~ "Eastern Europe and Central Asia",
                           e_regionpol_6C == 2 ~ "Latin America and the Caribbean",
                           e_regionpol_6C == 3 ~ "Middle East and North Africa",
                           e_regionpol_6C == 4 ~ "Sub-Saharan Africa",
                           e_regionpol_6C == 5 ~ "Western Europe and North America",
                           e_regionpol_6C == 6 ~ "Asia and Pacific",
                           T ~ NA_character_))

#------------- Figure 4: Protest against restrictions in different regions --------------------------------------------

# PLot regions
ggplot(acled_regions %>% drop_na(e_regionpol_6C), aes(x = reorder(factor(regio), restrictions_masks), y = restrictions_masks)) +
  geom_bar(stat = "identity", alpha  = .85) +
  ylab("Restrictions/masks (mean topic prevalence)") +
  xlab("") +
  coord_flip() +
  theme_bw()

ggsave(filename = "output/Figure_4.pdf",
       width = 7, height = 5, dpi = "retina")

#------------- Figure 5: Protest against restrictions in different regions --------------------------------------------

# aggregate at country-year level
acled_during <- acled_topics_vdem %>% 
  group_by(country, cowcode, country_text_id, year) %>% 
  summarise(across(c("health_care",
                     "restrictions_masks",
                     "vaccination",
                     "imprisonment_crime",
                     "mishandling_corruption",
                     "economic_consequences",
                     "business_restrictions",
                     "education",
                     "v2x_libdem"), 
                   ~ mean(., na.rm = T)),
            events = length(data_id)) %>%
  mutate(across(c("health_care",
                  "restrictions_masks",
                  "vaccination",
                  "imprisonment_crime",
                  "mishandling_corruption",
                  "economic_consequences",
                  "business_restrictions",
                  "education",
                  "v2x_libdem"), 
                ~ .x * log(events), 
                .names = "{.col}_weighted")) %>% 
  ungroup() %>% 
  drop_na(v2x_libdem) 

#-------------  Load and merge Oxford Covid data --------------------------------------------

# Add Coviod data
covid <- rio::import("data/owid-covid-data_update.csv") 

covid_agg <- covid %>% 
  mutate(cowcode = countrycode::countrycode(iso_code, "iso3c", "cown"), 
         year = as.numeric(str_sub(date, 1,4), .after = iso_code)) %>% 
  select(cowcode, year, everything()) %>% 
  drop_na(cowcode) %>% 
  group_by(cowcode, year) %>% 
  summarise(across(c("stringency_index"), 
                   ~ mean(., na.rm = T), 
                   .names = "{.col}_mean"),
            across(c("total_deaths_per_million",
                     "excess_mortality_cumulative_per_million",
                     "total_vaccinations_per_hundred"), 
                   ~ max(., na.rm = T), 
                   .names = "{.col}_max"),
            across(c("human_development_index",
                     "gdp_per_capita",
                     "total_vaccinations_per_hundred",
                     "hospital_beds_per_thousand",
                     "life_expectancy",
                     "population"), 
                   ~ first(na.omit(.)), 
                   .names = "{.col}_first")) %>% 
  ungroup()

# Record missing values
covid_agg[covid_agg == -Inf] <- NA
#covid_agg[covid_agg == NaN] <- NA

#merge with other data"
acled_vdem_corona <-  left_join(acled_during, covid_agg) %>% 
  drop_na(population_first) %>% 
  mutate(across(c("health_care",
                  "restrictions_masks",
                  "imprisonment_crime",
                  "vaccination",
                  "mishandling_corruption",
                  "economic_consequences",
                  "business_restrictions",
                  "education"), 
                ~ .x * (events/population_first), 
                .names = "{.col}_pop_weight")) %>% 
  mutate(tertiles_hdi = ntile(human_development_index_first, 3)) %>% 
  mutate(country_year = paste0(country, " (", year, ")"), .after = year)

# Scatterplots
min_label = .5
events_min = 1

# Development and restrictions
ggplot(acled_vdem_corona %>% filter(events > events_min),
       aes(x = human_development_index_first, y =  vaccination))  +
  #  geom_smooth(method = "loess", color = "black", alpha = .5) +
  geom_point(aes(size = events/population_first), fill = "lightgrey",
             colour="black", pch=21, alpha = .5) +
  geom_smooth(method = "lm", color= "black", alpha = 0.75) +
  xlab("Human Development Index") + ylab("restrictions/masks") +
  scale_size(trans = "identity", range = c(2,8), guide = "none") +
  ggrepel::geom_text_repel(data =  acled_vdem_corona %>%
                             filter(country_text_id %in% c("AUT", "CHE", "ESP",
                                                           "UGA", "USA", "COD",
                                                           "CHN", "RUS", "IND",
                                                           "VEN", "ITA", "TUR",
                                                           "BRA", "HTI") & year == 2021), 
                           aes(label = country), min.segment.length = 0) +
  theme_bw() +
  # ggpubr::stat_cor(method = "pearson", label.x = .5, label.y = .6) +
  theme(panel.grid.minor = element_blank())

ggsave(filename = "output/Figure_A3.pdf",
       width = 7, height = 5, dpi = "retina")


# Development and health care
ggplot(acled_vdem_corona %>% filter(events > events_min),
       aes(x = human_development_index_first, y =  health_care)) +
  geom_point(aes(size = events/population_first), fill = "lightgrey",
             colour="black", pch=21, alpha = .5) +
  geom_smooth(method = "lm", color= "black", alpha = 0.75) +
  xlab("Human Development Index") + ylab("health care topic") +
  scale_size(trans = "identity", range = c(2,8), guide = "none") +
  ggrepel::geom_text_repel(data =  acled_vdem_corona %>%
                             filter(country_text_id %in% c("AUT", "CHE", "ESP",
                                                           "UGA", "USA", "COD",
                                                           "CHN", "RUS", "IND",
                                                           "VEN", "ITA", "TUR",
                                                           "BRA", "HTI") & year == 2021), 
                           aes(label = country), min.segment.length = 0) +
  theme_bw() +
  #  ggpubr::stat_cor(method = "pearson", label.x = .5, label.y = .6) +
  theme(panel.grid.minor = element_blank())

ggsave(filename = "output/Figure_5.pdf",
       width = 7, height = 5, dpi = "retina")

#------------- Figure 6 and Tables A.2 - A.9: Regression analysis --------------------------------------------

acled_vdem_corona %>%
  select(health_care:education) %>%                   
  map( ~ lmer(.x ~ I(stringency_index_mean/10) + log(total_deaths_per_million_max) +
                I(human_development_index_first*100) + 
                hospital_beds_per_thousand_first + 
                log(population_first) + I(v2x_libdem*10) + log(events) + (1 | country) + (1 | year),
              data = acled_vdem_corona)) %>% 
  map( ~texreg::texreg(.x)) 

reg_results <- acled_vdem_corona %>%
  select(health_care:education) %>%                   
  map( ~ lmer(.x ~ I(stringency_index_mean/10) + log(total_deaths_per_million_max) +
                I(human_development_index_first*100) + 
                hospital_beds_per_thousand_first + 
                log(population_first) + I(v2x_libdem*10) + log(events) + (1 | country) + (1 | year),
              data = acled_vdem_corona)) %>% 
  map( ~ broom.mixed::tidy(.x, conf.int = T, conf.level = 0.95)) %>%  
  bind_rows(.id = "outcome")  %>% 
  filter(term!= "(Intercept)") %>% 
  filter(effect== "fixed") %>% 
  mutate(term = case_when(term == "I(stringency_index_mean/10)" ~ "Stringency index (avg.)",
                          term == "log(total_deaths_per_million_max)" ~ "Deaths per million (log)",
                          term == "I(human_development_index_first * 100)" ~ "Human Development Index",
                          term == "hospital_beds_per_thousand_first" ~ "Hospital beds per thousand",
                          term == "log(population_first)" ~ "Population (log)",
                          term == "I(v2x_libdem * 10)" ~ "Liberal democracy",
                          term == "log(events)" ~ "Protest events (log)",
                          T ~ "term")) %>% 
  mutate(outcome = case_when(outcome == "economic_consequences" ~ "economy",
                             outcome == "imprisonment_crime" ~ "imprisonment/crime",
                             outcome == "vaccination" ~ "vaccination",
                             outcome == "education" ~ "education",
                             outcome == "health_care" ~ "health care",
                             outcome == "mishandling_corruption" ~ "mismanagement",
                             outcome == "restrictions_masks" ~ "restrictions/masks",
                             outcome == "business_restrictions" ~ "business restrictions",
                             T ~ NA_character_),
         color = case_when(conf.low < 0 & conf.high < 0  ~ "sig",
                           conf.low > 0 & conf.high > 0  ~ "sig",
                           T ~ "nosig"),
         shape = case_when(conf.low < 0 & conf.high < 0  ~ "sig",
                           conf.low > 0 & conf.high > 0  ~ "sig",
                           T ~ "nosig"))

ggplot(reg_results, aes(x = term, y = estimate)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_point(aes(x = term, 
                 y = estimate, color = color)) + 
  geom_linerange(aes(x = term, 
                     ymin = conf.low,
                     ymax = conf.high, color = color),
                 lwd = 1) +
  scale_x_discrete("") +
  scale_color_manual(values = c("grey70", "darkblue")) +
  scale_y_continuous("Estimate") +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "none") +
  facet_wrap(~ outcome, nrow = 4)

ggsave("output/Figure_6.pdf", width = 8, height = 8, dpi = "retina")
